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We introduce a phase-space representation for qubits and spin models. The technique uses an SU(n) coherent 
state basis, and can equally be used for either static or dynamical simulations. We review previously known 
definitions and operator identities, and show how these can be used to define an off-diagonal, positive phase- 
space representation analogous to the positive P-function. As an illustration of the phase-space method, we use 
the example of the Ising model, which has exact solutions for the finite temperature canonical ensemble in two 
dimensions. We show how a canonical ensemble for an Ising model of arbitrary structure can be efficiently 
simulated using SU(2) or atomic coherent states. The technique utilizes a transformation from a canonical 
(imaginary-time) weighted simulation to an equivalent unweighted real-time simulation. The results are com- 
pared to the exactly soluble two-dimensional case. We note that Ising models in one, two or three dimensions are 
potentially achievable experimentally as a lattice-gas of ultra-cold atoms in optical lattices. The technique is not 
restricted to canonical ensembles or to Ising-like couplings. It is also able to be used for real-time evolution, and 
for systems whose time-evolution follows a master-equation describing decoherence and coupling to external 
reservoirs. The case of SU(n) phase-space is used to describe «-level systems. In general, the requirement that 
time-evolution is stochastic corresponds to a restriction to Hamiltonians and master-equations that are quadratic 
in the group generators or generalized spin operators. 



I. INTRODUCTION 

Qubits are a central concept in quantum information. How- 
ever, complexity issues mean that calculations with large 
numbers of qubits are nontrivial: the Hilbert space dimen- 
sion scales as 2 M for M qubits. A natural way to treat this 
type of complexity is to use a phase-space representation over 
an atomic coherent state basis. Coherent states, introduced 
by Schrodinger[l], have been used widely in quantum op- 
tics. Atomic coherent states - originally used for collections 
of two-level atoms[|2|] - are the natural solution for a quan- 
tum spin driven by an external driving force, like a magnetic 
field. They are also called SUY2) JHQ], spin, or more gener- 
ally SU(n) coherent states||5l, |7[] for arbitrary n-level systems. 
Since they are a continuous set, they satisfy differential iden- 
tities, which can have useful applications. 

In this paper, a phase-space representation of arbitrary den- 
sity matrices in terms of off-diagonal SU(n) coherent state 
projectors is introduced. This extends earlier P-function[6] 
and Q-function0 1^] approaches involving SU(2) and SU(n) 
projectors ! 1011 . The methods described here allow dynamical 
or static entanglement to be treated, and extend earlier phase- 
space approaches in quantum optics lfl2l [l3l [Til ITsl 11611. In 
particular, they include off-diagonal coherent-state projectors 
which lead to positive-definite diffusion, and hence to dynam- 
ical realisations as stochastic processes^ IH El. The re- 
sulting methods have applications to either time-evolution or 
canonical ensemble calculations of finite Hilbert space sys- 
tems with spin systems. More general applications in quantum 
information are also possible, owing to the simplicity with 
which large and/or decoherent spin systems can be treated. 

Other methods for treating finite Hilbert spaces like coupled 
spins include finite versions of the Wigner representation ll2lll. 
path-integral techniques l24ll and DMRG-based methods 11251 
I^IU While these are interesting and often very useful, they 
are not suited to exact, probabilistic simulations, because they 
either involve approximations, or else they do not use a posi- 



tive distribution function. When DMRG techniques are possi- 
ble - typically in one-dimensional ground-state calculations - 
they are very accurate and useful, but this method often can- 
not be used in many other physical examples involving finite 
temperatures, dissipation, dynamics or higher dimensions. 

Exact methods also exist - like the one and two dimensional 
Ising model at finite temperature - but these approaches are re- 
stricted to special cases. Our approach is to define a positive 
distribution function over a space of SU(n) coherent state am- 
plitudes. This is a much smaller dimension than the whole 
Hilbert space, scaling proportionally to the number of spins. 
We emphasize that the representation is not unique, and some 
care is needed in choosing the expansion to minimise sam- 
pling error. In general, the main restriction is the compactness 
or otherwise of the resulting phase-space distribution: if there 
are large distribution variances, this will increase sampling er- 
ror in a practical calculation. 

As an example to illustrate scaling behaviour in an ex- 
actly soluble case, the application of SU(2) or atomic coherent 
states to solving the two-dimensional Ising model is treated in 
detail. This application is simple yet instructive, and the re- 
sulting algorithm is novel and efficient. The Ising model l28ll 
is one of the oldest models in statistical mechanics, with many 
applications|29]. The model has the virtue of having a non- 
trivial exact solution in two-dimensions ll30l l3lll . It displays 
a critical-point phase-transition|32], which we use to test the 
phase-space method. We find excellent agreement with these 
exact results. 

The original use of the Ising model was a simple theory 
of ferromagnetism - in which atomic spins have either an 
'up' or 'down' orientation. It also finds applications to a va- 
riety of other physical problems, from the theory of lattice 
gases and binary alloys to spin glasses 0311 . percolation p"^ 
and other disordered systems. Modern ultra-cold atom exper- 
iments with optical lattices[35[] can test this model directly, at 
temperatures above quantum degeneracy where the lattice-gas 
model is applicable. In this case, the two states of each lattice 
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site correspond simply to the presence or absence of a single 
atom. At lower temperatures where coherences are important, 
Heisenberg-like models become applicable, and these will be 
treated elsewhere. 

There are numerous corresponding techniques for solving 
the Ising model. However, exact solutions are known only in 
special cases like the uniform one and two-dimensional lat- 
tices. More generally, the other techniques that are known 
rely on Monte-Carlo methods ll37l [38l [39ll . in which the space 
of all configurations is searched by random spin-flipping 
algorithms 1 40]. The method demonstrated here is quite dif- 
ferent to traditional approaches. 

The SU(n) phase-space approach can also be readily used 
for other models of interacting spins, to real time evolution 
and to dynamical couplings to reservoirs, where no exact so- 
lution is known. While these applications will be treated else- 
where, we note that the main restrictions are that the Hamil- 
tonian or master equation should be at most quadratic in the 
SU(n) operators, which is the typical case for coupled spin 
systems. It is intriguing to note that these types of prob- 
lems are also regarded as potentially soluble for future gen- 
erations of quantum computers. The methods proposed here 
have the advantage that they can be implemented on digital 
computers. Thus, they complement the quantum computing 
approach, and indeed can be used to simulate quantum logic 
gates in the presence of decoherence. The main limiting is 
sampling error, which typically grows with simulation time. 



II. SU(2) COHERENT STATES 



2S spin 1/2 quantum systems or qubits, each with 2-levels 
and equivalent couplings. These composite systems in general 
have 25+1 distinct energy levels, and there is a unique lowest 
eigenstate of S z , denoted |0). 

The standard definition of SU(2) coherent statesiHIIl] is that 
they are the states generated from |0) by the raising operator, 
so that, for a spin-5 basis, 



l«) (2) = ' 



7 10) 



(2.4) 



It is convenient here to also consider an un-normalized ver- 
sion of this atomic coherent state, which we define as 



\\f)W = [y°] N eV 1S+ ^\0) 



(2.5) 



For simplicity in obtaining identities, it is useful to have just 
one complex parameter, as in the standard definition. Our 
choice is to define 



y/ 1 = exp(z/2) 
y/° = exp(-z/2), 



(2.6) 



where z = r + i(j) = In a is a complex parameter. With this 
choice, the SU(2) coherent states are parametrized over a 
one-dimensional complex manifold, or a two-dimensional real 
manifold. We will represent this parametrization as \\z), where 



We start with the well-known SU(2) case, which corre- 
sponds to a spin-J physical system or more generally, a 
collection of physically equivalent two-level systems. The 
SU(2) coherent states or atomic coherent states are defined 
for states generated with angular momentum raising and low- 
ering operators(3|,|4j]. These are physically important in many 
systems, ranging from groups of two-level atoms to nuclear 
spins, as well as superconducting qubits and other systems 
with an SU(2) symmetry. 

The relevant spin operators S have commutators defined so 
that: 

[Sjj]=ie ljk s k . (2.1) 

Here e i; <. = ± 1 depending on whether the indices are in cyclic 
or anti-cyclic order, and one conventionally writes §x,y,z to de- 
note 5i ,2,3- It is useful to also define the raising and lowering 
operators which act on an eigenstate of S z to increase (de- 
crease) the eigenvalue. These are defined as: 

S ± =S x ±iS y . (2.2) 

We consider a subsystem with a definite value of 

S 2 x +S 2 + S 2 = S 2 = S(S+l). (2.3) 

Physically, these may be obtained either directly as an atom 
or molecule of spin S, or equivalently from a grouping ofN> 



\\z) = e- s 'e s+eZ |0) . (2.7) 

For visualization purposes, one may project the atomic co- 
herent state phase-space onto a spherical surface, called the 
Bloch sphere. In this case, it is usual to normalize the state, 
and to define 

|0,0> (2) = |e^tan0/2) (23 . (2.8) 

This Bloch-sphere mapping therefore involves the transforma- 
tion of 

a = entail 0/2 = e' . (2.9) 



1. Two-level case 

As an illustration of the simplest case possible, where 
S = 1/2, we consider a two-level Hilbert space having quan- 
tum states labelled |0) and |1). This corresponds to a single 
qubit in quantum information terminology. An atomic coher- 
ent state or SU(2) coherent state is then just an arbitrary pure 
qubit state: 

||V) (2) = ^°|o) + r 1 |i> 

= e- z l 2 \Q)+e z l 2 \\) . (2.10) 
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This shows the utility of this parametrization: it displays a 
symmetry between up and down states, which simply corre- 
sponds to changing the sign of z. In a useful vector represen- 
tation, one can write this in an explicit form as 



\w 



,(2)_ 



(2.11) 



In this notation, the state |1) corresponds to spin projection 
m = 1 /2. Similarly, the second entry or state |0) corresponds 
to spin projection m = — 1/2. On the Bloch sphere, this cor- 
responds to 



|0,0> 



(2) 



v°\o) + v l \i) 



0|2 



V 



= cos — e 
2 



sin — e 
2 



i<j,/2 



(2.12) 



1) (2.13) 



A. Lattice atomic coherent states 



which are generated using operators with an SU(n) operator 
algebra. 

The SU(n) group is the group of n x n unitary matrices with 
unit determinant, and so provides the most general way to treat 
the transformations of an n-level quantum system. There- 
fore, SU(n) coherent states provide a useful basis set for gen- 
eral multi-level quantum systems like atoms or spins. In the 
following section we review results for the SU(n) coherent 
states. We also consider the important case of outer products 
of SU(n) coherent states, which are needed for treating lat- 
tices. 

In the simplest case n corresponds to the number of distinct 
quantum states or levels involved. More generally, n simply 
labels a symmetry group which can have a larger dimensional 
representation, just as in the SU(2) case. 

These states are useful in treating, for example, an assem- 
bly of n coupled Bose-Einstein condensates, n-level atoms, 
or photon states with 0,1 ...n — 1 photons per mode. The 
SU(n) algebra is generated by the n 2 — 1 independent oper- 
ators which satisfy the commutation relations[41] 



For M distinct spins, particles, or lattice sites, where one 
may wish to address or couple to them individually, one must 
have M distinct spin operators. As noted above, each of these 
can describe N physical qubits. 

The corresponding outer-product SU(2) coherent state is 
then: 



a 



(2.M) 



M 

n 



1+ cu 



|0) 



(2.14) 



For N = I, the two-level or qubit case, we note that with z 
Z = (zi-.-Zm), our un-normalized definition becomes; 



(3.1) 



together with the constraint that = 1. The SU(n) co- 

herent states can also be written in the following convenient 
form, using an un-normalized notation in analogy to Eq. (12.5b . 

as: 



\\yt)M = [ W } N e^>°^°/V \0) 



(3.2) 



We can use a collection of N equivalent n-level quantum 
systems with states ■ for ji = 0, ...,« — 1, and j = 1, . . .N , 
to indicate the essential features of this approach. In this case 
the SU(n) operator algebra representation is provided by: 



(2,M) 



J m=l 



,-z,„/2 



|0)„ 



n/2\ 



(2.15) 
(2.16) 



In this notation, the inner product is 

M 

-,M 



z||||z')=2 M n cosh ([4 + 4]/2) , (2.17) 



= i 



and we can therefore introduce a normalized state denoted |z), 
where 



z = 



M 

n 



=1 v / 2cosh(r,„) 



(2.18) 



III. SU(n) COHERENT STATES 

In cases where SU(2) symmetry does not hold, the SU(2) 
coherent states can be generalized to SU(n) coherent states 



(3.3) 



For this case of N equivalent n-level atomic or spin states, 
one can then define an SU(n) coherent state directly in terms 
of the original Bloch basis \k) as: 



» 



N 

n 

.7=1 



(3.4) 



The corresponding normalized state is then: 
1 



X;VlM>; 

M=0 



(3.5) 



In the normalized case it is common to take the first coef- 
ficient to be unity, so that y/° = 1, although other choices are 
possible. In general there are n — 1 independent complex am- 
plitudes of physical significance, since the overall phase and 
amplitude of a wave-function is physically irrelevant. 



4 



A. Lattice SU(n) coherent states 

Lattice coherent states we introduced in a pioneering work 
of Shastri et aljlt], to study the Heisenberg model of interact- 
ing spins. In our notation, for SU(n) coherent states defined at 
multiple sites on a lattice labelled m, we introduce: 



[n.M) 



= lllV»; 



(3.6) 



or, in a matrix notation analogous to the two-level case - ex- 
cept with n levels per site - 



(n.M) 



n-\ 
m 

-n-2 



Wm 



(3.7) 



These multiple SU(n) coherent states have the following inner 
products: 



,\ (n,M) M r—^ 
m— 1 



(3.8) 



One can also introduce normalized SU(n) lattice coherent 
states, where the normalization uses the distance measure 



N = I, this can be simplified further, as one obtains from the 
z-parameter mapping that 



~ _ r 2n de 

Jo 2k 



i«/2|l) + < .-»/2|o)l \ e -°/ 2 (l\ +e ' e / 2 (0\ 



271 dd „ w , 

x- 10/0 
o 2% 

271 d 9 , w , 

— iB)(i6 
o % 



(4.2) 



In the more general SU(n) case, one finds thatEdl: 



i = / 5 (ivi 2 - i)iy> (VM 2n v. (4.3) 

An even simpler resolution of the identity operator (for = 1) 
is easily obtained with a multiple phase integration: 



1 



1>(1| + |2) (2| + ... + |„>(»| 

In rln d n-\Q ^ 



2n 







o {2nf- 1 

271 d"- l e 

(27T 



,n-l 



(4.4) 



Just as in the two-level case, the first phase integral is omitted 
here (ie, 9q = 0), since this term is always orthogonal to the 
others, due to the remaining phase-integrals. 



Hence: 



(n.M) 



M 

n 



1 Wm 



(3.9) 



(3.10) 



These kinds of states can be thought of as generalizations 
of the harmonic-oscillator coherent states, in the sense that 
with the usual harmonic-oscillator coherent states there are 
prescribed relationships between the coefficients. In the SU(n) 
case there is no fixed relationship between coefficients, but 
there is a fixed upper bound to the quantum number. 



IV. COMPLETENESS AND IDENTITIES 

A. Completeness 

The spin coherent states form an over-complete basis. In 
the SU(2) case with spin-S, the resolution of the identity is 
well-known(4], and is given by 



1 = {2J+V 



dQ. 
An 



(4.1) 



where dQ. = dcos 6d(j) is the usual integration measure for the 
solid angle in spherical coordinates. In the spin-half case with 



B. SU(n) operator identities 

We wish to obtain differential identities that involve the set 
of operators that can act on the spin coherent states. These can 
all be regarded as extensions of the very simple differential 
identities that exist for the SU(n) coherent states. From Eq 
(13. 4t , one can directly prove that: 



R. 



'n.M) 



Wn 



dy/n, 



(n.M) 



(4.5) 



We now specialize to the two-level case where 'raising' and 
'lowering' operators are conventionally defined in physics as 
the matrices: 



(4.6) 



These have a direct relationship with the R operators, since for 
SU(2) symmetry with S = 1/2, one has: /? 01 = <J- and = 
(7+. In addition, a x,y ' z are the Pauli spin operators defined as: 



"0 1" 




0" 





; a = 


1 



"0 1" 




-z" 




' 1 


1 


, d y = 


i 


, s z = 


-1 



(4.7) 



Here as well, there is a correspondence with SU(2) genera- 
tors, since 



1 



(a x ±ia x ), 



(4.8) 
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and 



(4.9) 



Identities can either be obtained from these correspondences, 
or from direct differentiation, since: 



(4.10) 



d 


\ ^ ' 


1 


' e'l 2 ' 






~ 2 





Hence, in operator language: 

§- z \\z) = ~c*\\z). 
On taking the hermitian transpose: 

4i;(z\\ = Uz\\o z ■ 
dz* 2 x " 

With a little algebra, one can also show that 



(4.11) 



(4.12) 



a positive-definite diffusion and hence to stochastic equations 
that can be numerically simulated. 

Instead, we will focus on the SU(2) and SU(n) cases anal- 
ogous to the positive P representation !!?! [Till . This approach 
includes off-diagonal projection operators in the expansion 
of the density matrix, and give rise to a phase-space dimen- 
sion which is at least twice that of the classical phase-space. 
The result is a complete, positive representation that gener- 
ates positive-definite Fokker-Planck equations. T his g eneral- 
izes related work in quantum and atom optics|22, 23], which 
uses similar procedures. 



A. SU(2) phase-space expansions 

We now illustrate these ideas with reference to the simplest 
SU(2) or qubit case, using the reduced z— parametrization. 
If the density matrix is separable, one can use a representa- 
tion in terms of a positive probability over the SU(2) diagonal 
coherent-state projectors: 



1 d_ 

2 Tz 



oHz) 



(4.13) 



C. Equivalent Identities 



( z )|z) ( Z |rfz. 



(5.1) 



It is always possible to define a positive representation like the 
Husimi Q-function, which is: 



Here the functions differentiated are all analytic functions, 
either of z or of z* ■ This means that we can always use 
Cauchy's equivalence of differentiations in real and imaginary 
directions, i.e., 



dz 
_d_ 

dz* 



z = 



— II ) = — 

dr llZ> d(j) 



dtp 



(4.14) 



This freedom, which also applies in the SU(n) case, allows 
one to derive a variety of different equivalent equations for a 
given operator evolution equation. 



V. SU(N) PHASE-SPACE 

Just as with the harmonic-oscillator coherent states, it is 
possible to define a variety of operator representations using 
the SU(n) coherent states. A number of these have been ex- 
tensively studied, including representations analogous to the 
Wlfl2tl, Olfll, P0], and +Plfl7|jl8tl representations. Spin 
versions of the Q-representationQ, P-representationQ and 
Wigner representations !^ have been introduced previously. 
These essentially are defined on classical phase-spaces, in the 
sense that the phase-space dimension is the same as that of the 
generators of the coherent state. 

However, as in the case of the harmonic oscillator, these do 
not generally allow time-evolution equations with a stochastic 
(positive) propagator. The difficulty here is that in general, 
these types of phase-space representation do not give rise to 



Q 



(2), 



MP * 



(5.2) 



However, these two methods will not generally give a 
positive-definite diffusion in the time-evolution equations for 
the distribution, except in special cases. In order to achieve 
this, we must introduce off-diagonal coherent state projectors, 
resulting in an expansion of form: 



(5.3) 



Here we define A = (z,z'), so that dX = d 2N zd 2N z', and we 
have introduced a general kernel operator with an arbitrary 
weight coefficient w: 



Al 2) (A) 



Xl 2) (z,z') = ||z) (z'\\e^ 



(5.4) 



With the simplest choice of w = 0, we obtain an expansion 
in terms of un-normalized projectors, which from Eq (12.17b 
leads to the result that 



p = J ( z ,z') A<, 2) (z,z') d 2N zd 2N z' , (5.5) 
with a trace given by 

Tr(A< 2 )(z,z')) = (z'\\z) 

= n[2cosh([z}+4]/2)] 

7=1 

= A(R), (5.6) 
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where we have introduced the kernel trace A (R) as a function 
of the combined variable R = [z* + z'] /2. 

There are many other choices of weights and phase-space 
expansions. One choice is to define the weight w(z,z') = 
— ln(z' ||z). This choice ensures that the kernel has a unit 
trace, giving results analogous to the positive-P approach. In 
this case: 



Xt 2) (A) 



:(2) 



(z,z') = 



(5.7) 



More generally, either using Ao as a dynamical variable, 
or other choices of weight function are necessary, in order 
to eliminate boundary terms which can arise in dynamical 
equations tl2(X 14311 . 



P = Jp { " ] (^_)\f){f\d^_. (5.11) 

A positive Q-function like phase-space representation always 
exists, with: 



(5.12) 



Just as in the SU(2) case, neither of these phase-space meth- 
ods will usually result in positive-definite stochastic evolution, 
either for canonical ensembles or for dynamical evolution. To 
overcome this limitation, a positive representation using off- 
diagonal projectors must be introduced: 



B. Entanglement and Bell states 

We note here that there is a fundamental contrast between 
this approach and the diagonal P-representation approach 
originally due to Sudarshan and Glauberjl4|, and later ex- 
tended to SU(2) coherent states[4]. The basis set of the diago- 
nal P-representation is separable: it therefore cannot represent 
entanglement, except as a limit of a generalized function. 

By comparison, the present approach includes terms that 
are fundamentally inseparable, and therefore can represent 
states like Bell states. To see this, consider the Bell state de- 
fined as: 



\r) — =[|o,i>- 



V2 

i 

1 

V2 



1,0)] 

r 



(5.8) 



where: 



1 H 



1 o 
1 

1 

1 



(5.9) 



The corresponding density matrix is: 



(5.10) 



This has the form of a positive distribution over the off- 
diagonal coherent state basis terms, as required. 



C. SU(n) phase-space expansions 

We now consider the most general SU(n) case. It is well- 
known[7] that one can define a diagonal phase-space repre- 
sentation analogous to the Glauber P-function: 



p= J pW(X)A w W (X)dX. (5.13) 
Here we define X = (Ao,~y?,~0), so that dX = d 2 ( d+ ^X = 



d 2 Xod 2M "y}d 2M " (j) where d = 2Mn, together with a general 
kernel operator with weight coefficient w: 



a«(A)=aW(v^ 



(n) 



, x (»> M ) , 



(h,M) 



(5.14) 

This reduces to the diagonal case when ~y? = . From Eq 
(|2T7j, the simplest choice of Ao = w — leads to the result 
that: 



Tr(A«(T?,? 







(«,M) 



(n,M) 



n [c-v^» 



Another choice is to define the weight 



w 



= -ln $ 



(5.15) 



(5.16) 



so that the kernel has a unit trace, giving results analogous 
to the positive-P approach. However, unless there is damping, 
this choice by itself can lead to instabilities and boundary term 
errors JH. 

If Ao 0, it can be used as another dynamical variable, giv- 
ing stabilized weighted trajectories as in the stochastic gauge 
method[20]. More general weight choices are also possible. 
The use of different weights changes the form of the resulting 
dynamical equations, thereby giving rise to useful techniques 
which can be utilized to optimize and solve these equations. 
An example will be given in the next section. 



VI. DYNAMICAL CALCULATIONS 

The calculation of observables and correlations in real or 
imaginary time (for thermal equilibrium) is the main purpose 
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of this phase-space method. The advantage of the approach is 
that it is a general-purpose method. The identities and trans- 
formations involved do not depend on detailed properties of 
the Hamiltonian, apart from the requirement that it must be 
able to be expressed using the group generators. 

Provided this requirement is satisfied, the calculations in- 
volved are not specific to a given model. However, some cau- 
tion is necessary. The probability distributions obtained can 
have a variety of widths in phase-space, which means there is 
a large range of potential sampling errors possible. This is not 
uniquely specified by the Hamiltonian. As the SU(n) basis set 
is not orthogonal, the phase-space distribution is therefore not 
unique, and depends on the precise identities and algorithms 
chosen. Since the underlying coherent states factorize on a lat- 
tice, one may expect that increasing correlations and entangle- 
ment between lattice sites will require an increased 'footprint' 
of the distribution, and hence an increased sampling error. 



A. General evolution problems 



can then be transformed into the stochastic equations, which 
in Ito calculus are generically of the form: 

dXo/dt = U + g^ li --g ll gf 1 

dX^/dt = A m +B mv (Cm-^)- ( 6 - 5 ) 

Here the weight term U and the drift vector A are deter- 
mined by the form of the original Liouville equation. The 
drift gauges appear as the arbitrarily functions g, and diffu- 
sion gauges appear as the freedom that exists in choosing the 
noise matrix B. The noise terms £ are Gaussian white noises, 
with correlations: 

(C M (t)Cv(0> = fyv 5 ('-0- (6.6) 

Equations (l6.5H6.6t can be used to solve a large class 
of quantum dynamical and thermal-equilibrium problems in 
coherent-state representations. In practice, the numerical im- 
plementation of these equations can be simplified by use of 
automatic code-generators Pfl l45ll . 



To illustrate the procedure, the required dynamical evolu- 
tion is first written as a Liouville equation for the density op- 
erator. This may or may not be unitary, and does not have 
to be trace-preserving, as long as it is linear in p, and can be 
written using a polynomial in the group generators: 



dp(t)/dt = L[p(r)]. 



(6.1) 



To solve this with phase-space methods, we first expand the 
density operator over the SU(n) operator basis A^ M \X), 
where X is the set of all complex coherent amplitudes: 



p{t) = j P(X \t)A{X)d~X 



(6.2) 



This defines a d + 1 -dimensional complex phase-space, where 
d = 2Mn as before, with a dynamical weight variable Xq if 
necessary. The SU(n) differential identities allow us to write 
the Liouville operator equation as 

dp(t)/dt = f P(t,t)£> A A(X)d 2{d+l) t, (6.3) 



where Jz^ is a linear differential operator. Due to the non- 
uniqueness of the identities, this can include arbitrary stochas- 
tic gauge functions. Provided there are no derivatives higher 
than second order, this equation can finally be transformed 
into a positive-definite, weighted Fokker-Planck equation for 
P. It is essential that the gauges are chosen to eliminate 
any boundary terms that may otherwise arise from the partial 
integration SEl- 



U-—A 1 92 D 
dXp 2 dX^dX v ^ v 



P(X,t). 



(6.4) 

Here we use a summation convention where pL — 1 , d. In- 
troducing a matrix square root B, where D^ v = B^pByp, this 



B. Operator identities: SU(2) case 

To use this approach, one must obtain differential identities 
for the group generators. We start with the SU(2) case. Here 
we will omit the superscript (2) indicating an SU(2) kernel, 
when there is no ambiguity. 

With the simplest constant weight choice we will use here 
of w — 0, the only differential identities needed are obtained 
directly from Eq ( 14.1 U and Eq ( 14.12b i.e., 



S Z A , 



3— Aq = A S Z . 
dz* 



(6.7) 



Other useful differential identities in more general cases are 

A w , 



d A 

t-A w = 
dz 



dw 

S ~ + lTz 



dz* 



A w 



S z + 



Hence, for example, one can write: 

c<X <5 dw 

o A w = — r 

dz dz 

A W S Z = 



dz'* 



A w , 



(6.8) 



(6.9) 



dz 1 * dz'* 



C. Operator identities: SU(n) case 

We wish to obtain similar differential identities for the 
SU(n) coherent state kernels. These are: 



V, 



1\\V 



dWm 



ftllv ,, dw 



8 



d<j>, 



— A W - A W 



fx * ' 



dy/(pn, 



(6.10) 



Since each occurrence of a group generator R gives rise 
to a differential term, the requirement that time-evolution is 
stochastic corresponds to a restriction to Hamiltonians and 
master equations that are quadratic in the group generators 
or generalized spin operators. 



D. Observables 

We illustrate how to calculate observables by reference to 
the the spin-half system , where the main observable of inter- 
est is the magnetization at site i, given by: 



Trip) 



(6.11) 



Defining the normalization as Z = Tr(p), with a measure 
dX = d 2N zd 2N i! ', one obtains that the uniform weight expan- 
sion case has the normalization 



Z = j P(z,z')A(R)rfA , 
= <A(R)) P . 



(6.12) 



Noting that 



/V 



Tr(a?A ) = tanh (/?,•) f~[ [2cosh(7? ; )] , (6.13) 

we can introduce a c-number equivalent magnetization vari- 
able nij = tanh (Rj). The mean magnetization is then written 

as 



P (z,z') tanh (/?,•) A (R)<a 



= -(tanh(^)A(R)>i>- (6.14) 

Similarly, the correlation function between two different 
sites is 



af a]) = - (tanh (J?,) tanh (Rj) A(R)) l 



(6.15) 



E. Phase-independent case 

In the case where the Hamiltonian is only a function of a z 's 
- as in the Ising model, considered in the next section - a 
much simpler expansion of the density operator can be used. 
While this is less general, it provides an alternative way to 
derive the results in the next section. 

This simplified expansion is: 



P(R)A,(R)dR, 



(6.16) 



where A z (R) is obtained on phase-averaging over the com- 
plete kernel, with the result that: 



M 



a z (r) = n^P' A, ^/ : 

7=1 



M ~ iRj(T Z 



m 

7=1 n=0 
M 



Y\2 (cosh{Rj) + djsmh{Rj)) . (6.17) 

7=1 



The operator correspondence 



07A- (R) = —A, (R) 

7 ^ ' dRj ZK ' 



(6.18) 



then holds. 

In the following section, we will focus on using the full 
coherent state identities, as these are more generally applica- 
ble. However, we note that for those primarily interested in 
the Ising model, our results can also be readily obtained using 
this reduced expansion. 



VII. THE ISING MODEL 

As an instructive example, we show that a lattice of SU(2) 
coherent states can be used to solve for the partition function 
of the Ising model of interacting spins. This is the simplest 
nontrivial case where one obtains an exactly soluble phase- 
transition in a spin model in two dimensions. As well hav- 
ing a wide applicability, it does illustrate many of the funda- 
mental scaling issues that occur in using phase-space methods 
to solve coupled spin models. Similar features also occur in 
more complex quantum spin models, which will be treated in 
greater detail elsewhere. 

Although we focus here on the simplest case possible where 
S = 1/2 at each site, we note that the basic ideas also hold 
for more general coupled spin-S 1 spin systems, or interacting 
atoms described by the most general SU(n) coherent states. 
However, in this example we make use of some identities and 
simplifying features that are unique to the spin-half case. 

The most general form of this model - in a summation con- 
vention which sums repeated indices - has the Hamiltonian 



H = -£/2;Sf--£Ji/Sf3j 



(7.1) 



We will assume here that the coupling term /y is symmet- 
ric, with > 0, which corresponds to attractive interactions 

between spins. Since (ofj = 1> self-interactions have no ef- 
fect apart from shifting the energy origin, and therefore it is 
common to set = for simplicity. Different choices of Jij 
will generate different types of Ising model, that can have any 
dimensionality, shape, or distribution of interaction strengths. 
The choice of /y = J for all nearest neighbours corresponds to 
the standard Onsager model lBOIl . The interaction terms will be 
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called links, since they typically join neighbouring spin sites 
or nodes on a lattice. The factor of half in the Hamiltonian ac- 
counts for the fact that all links are counted twice in the double 
summation. 

The density matrix, which gives information about the spin 
distribution in thermal equilibrium, is 



p = exp ( -j5H 



(7.2) 



One often wishes to calculate the total partition function Z, 
where 



Z = Tr(p). 



(7.3) 



If all the terms Jjj are either equal to each other or zero, then 
the interactions are uniquely characterized by a graph showing 
which nodes are linked by a nonzero interaction. Hence, there 
is a close relationship between the Ising model, and mathe- 
matical problems that count paths on a lattice. Once the total 
number of ways of constructing links with a given energy is 
known, the partition function can be easily obtained. Since 
there are 2 M distinct spin configurations, it is exponentially 
difficult to evaluate this directly, unless special types of sym- 
metry occur which can sometimes lead to exact solutions. Ex- 
amples are the case of the one and two dimensional regular 
lattice with uniform nearest-neighbour interactions, and the 
simplex with all node-pairs linked equally. 

More generally, one must use probabilistic methods to sam- 
ple the spin configurations. The standard techniques involve 
Monte Carlo or Metropolis techniques in which spins are 
flipped randomly, in order to obtain an ensemble of spin con- 
figurations at a fixed temperature. There is a long history to 
these methods, which can give excellent results. However, 
while much more efficient than direct configuration counting, 
these methods are still computationally intensive. This means 
that there are often strong limits to either the size of the lattice 
or to the accuracy, which is limited by the sampling error. Re- 
cent improvements in these standard techniques involve flip- 
ping clusters of spins, which is more effective at the critical 
temperature where the correlation lengths are large. 

We consider a different approach to this calculation using 
a differential equation method, that uses the atomic coherent 
state basis with a continuous parameter, rather than discrete 
spin configurations. The density operator satisfies the follow- 
ing equation ll42ll : 



dp 



The initial condition at high temperature is just 

p(o) = i 



(7.4) 



(7.5) 



A. Fokker-Planck Equation 



From the two-level completeness identity, Eq ( 14.21 i. one can 



write 



M 



p(zy,o)=n 



-5(9,-9!) 8^)8^ 



(7.7) 



This involves a single unique r value, rj = 0, and a random 
phase. This is transformed using operator identities into the 
resulting Fokker-Planck equation is transformed to a stochas- 
tic differential equation that can be sampled. We can choose 
equations in which the initially random phase is invariant. 
This leads to a stochastic equation in rj in which the initial 
state is given exactly, without sampling error. This technique 
can also be written as a type of path-integral. 

To illustrate the idea, we start with the simplest unweighted 
kernel, as previously: 



(7.8) 



p(/3) = J P(X,P)A(X)dX 

= Jp(z,z\p)\\z)(z'\\dX 

We see from this that 
dp(P) 1 



P(k,P) A(X),H 



dX. 



dp 2. 
Introducing the mean interaction strength per spin, 
1 



J 



2N 



(7.9) 



(7.10) 



we then rewrite the Hamiltonian in a form that allows us to 
obtain positive-definite diffusion terms, 



1 



H = -hiGf--J, 
= H'+MJ 



-NJ 



(7.11) 



The constant term has no effect on observable quantities, 
and will be neglected in the following calculations. In other 
words, we will calculate 



exp(-pH') =pe 



JMJ 



(7.12) 



which differs from the p defined above only by an overall 
normalization factor. Inserting the relevant identities, the two 
different operator orderings give 



l ~H>A 
2 



A, 



1 ? 

h i d i +-J ij (d i + d j Y 



(7.13) 



and: 



Next, the partition function p is expanded using an SU(2) 
coherent state projector basis, so that 



p(j8)= P(z,z',p)A{z,z')dX 



(7.6) 



AH 

2 



l,,a; ■ -J, J [a; ■ a, 



hi d' i + l -J ij (d' i + d'jf 



A. (7.14) 
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Here we have used the definitions <9, = <9/<9r,- and <?', = 
d/dr'i. We now introduce an extended vector notation with 
indices /I = 1, . . . ,2N, so that r^+j = r'j and dp = d/dr^ , 
with coupling constants J^ v , defined so that Ji+N.j+N = Jij, 
and hi + N — hi . 

Next, on integrating by parts, and equating coefficients of 
A, one obtains the following Fokker-Planck equation, with ex- 
plicitly positive definite diffusion terms: 



dP_ 
d/3 



1 - r 



P. (7.15) 



B. Stochastic Equation 

To obtain an equivalent stochastic equation, we must first 
write the Fokker-Planck equation in the form: 



9P rl 



Afi + —Dfi V d v 



(7.16) 



A suitable factorized diffusion matrix form is readily found 
by expanding the diffusion matrix D^v as a sum over distinct 
terms for each non-vanishing link, that is: 



2N 

fi,V 



[. . . lju, . . . , ly, • • •] • 



(7.17) 



It is immediate that D can be factorized in the form: 



2N 



£=^/ Atv B^ v) B^ v » r , 



(7.18) 



where B^ v ) is a 2N dimensional vector with two non- 
vanishing entries at and v respectively, i.e., 



g(Mv) 



The corresponding stochastic equations are then: 

dr 2N 



^JL - A + V R (M ' V V , , 
= V+E(^v + Cv M ) • 



where the independent real stochastic noises are corre- 
lated as 

5(j3-/3')- (7.21) 

These equations have the feature that they involve noise 
terms that are automatically correlated between pairs of spins 
linked by an interaction term, Jjj. The initial random phase is 
not changed by the interactions, and only the magnetization 
- which depends on rj - changes randomly in time. Spins 
that are linked tend to change together, as they experience a 
correlated noise term. 

Only the sum of r 7 + r' ; = 2Rj is relevant to the observed 
spin orientation. Defining 

Wtj W) = \j* («y (£') + Q + n, j+ n (p'))dp , (7.22) 

the resulting noise terms have a variance proportional to the 
inverse temperature: 



C. Partition function 



(7.23) 



The solution at inverse temperature j3 is: 



(7.24) 



where W,.+ (J3 ) = W tj (J3 ) + W n (j3 ). The resulting partition 
function is simply obtained on averaging over all the stochas- 
tic trajectories, so that: 

Z(j3) = (A(R{P)))e-P Nj 

= (^\[2 C o S h(Ri(P))]y~ PNf - (7-25) 

This gives an explicit solution for the partition function as 
an expectation value over the random processes R(j3). We 
note that while one may try to evaluate the partition function 
by simply averaging over many stochastic trajectories, this is 
far from being an efficient procedure. The problem is that 
the weights A(R(j3)) grow exponentially large for large val- 

(7.19) 

ues of \Rj\, which results in a large dispersion of trajectory 
weights, and therefore extremely large sampling errors. This 
naive method is not practical. A much more efficient proce- 
dure will be given in the next section. 

We notice at this stage, however, an interesting feature of 
these results. This is that the noise terms act only to couple ad- 
jacent sites together. Thus, an understanding of the renormal- 
ization behaviour of this problem can be realized by grouping 
spins together into clusters, in which case the residual noise 
(7 20) fr° m cluster interactions scales proportionate to the surface 
area of the cluster, rather than from the total volume. 
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1. Example: 2-site problem 



VIII. COMPUTATIONAL STRATEGIES 



As an example of the simplest nontrivial case with a uni- 
form external field (i.e., /i, = h) the two-node partition func- 
tion has only one link, so 



H = -h (of 



Ja\a : 



(7.26) 



There are four distinct states with interaction energies of ±7. 
Taking the trace, one can directly check from expanding over 
the four-dimensional configuration space, that 



There are several possible strategies for calculating the par- 
tition function while taking account of the final weight. For 
the Ising model, a direct solution to the original stochastic 
equation is inefficient for large M, as almost all trajectories 
will have an exponentially small weight compared to a very 
small number of optimal trajectories. We will demonstrate a 
strategy for making use of the fact that we now have a solu- 
tion to the stochastic equations in closed form, which allows 
the problem to be re-sampled in a more efficient way. 



Z 2 = Tr 



-pH 



= e P(2h+J) +2e -jp +e P(-2h+J) ^ (7 2?) 

For h = 0, the two-site correlation can be calculated imme- 
diately to be 



tanh (J37) . 



(7.28) 



We now wish to demonstrate how identical results are 
obtainable from the raw stochastic equations. Introducing 
W + (]8) = Wn (J3) + W 2 \ (j8), with (W +2 (J3) > = j8/, one finds 
that the two SU(2) coherent state amplitudes are always equal 
to each other: 



/?i(/3)=* 2 (/3) = /ij3 + W + a3) 



(7.29) 



Hence with J = 7/2, the partition function calculated from the 
stochastic equations is 



ZQ3) = {^\[2cosh(R,(l5))}j e 

2 



-fiNj 



-pj 



(7.30) 



Now, for a Gaussian process, 



,±2*03) 



(7.31) 



so the final result for the partition function is 

Z(J3) =e P(2h+J) +2e -jP +e P(-2h+J) , (7J2) 

Similarly, for the correlation function in the limit of h = 0: 



< Sf6 ^ = |7jy< tanh2 W A ( R ))/ 



4e -pJ 



(sinh 2 (*)>, 



Z(P) y "' /p 
= tanh(j3/). 



(7.33) 



This agrees with the result from the direct calculation. 



:exp[±2/z/3+2(W 2 (j3)) /) ] = exp[±2/zj3 +27)3] , 



A. Optimized stochastic methods 

One way to solve this problem is to use weighted kernels or 
gauge equations, combined with a strategy for breeding tra- 
jectories of largest weight, which is essentially the diffusion 
Monte-Carlo approach II 3 711 . Another approach is to use the 
Metropolis method l40ll . in which the link noise Wy is repeat- 
edly randomized, based on the final weight it generates, with 
some choices being accepted and some being rejected. 

A third way is to define a new stochastic equation whose 
solution gives the link noise distribution, without any addi- 
tional weight. To see this more clearly, suppose we write the 
final partition function as a multi-component integral over the 
link noises W, including the Gaussian weight factor used to 
generate the noises Wif. 

Z(/3) = y...y</Wexp(-V(W,j3)) , (8.1) 

where we have ignored all irrelevant normalization terms, and 
introduced a potential that already includes the weight factor: 

v ( w '£) - L j^ w u - D ncosh + Wj (0)) • 

(8.2) 

The first term is the most important at high temperatures. It 
tends to keep all link noises small, so that the magnetization 
is nearly zero. The second term is increasingly important at 
large j3, as it gives an increasing weight to terms with large 
correlated noises Wij, in which all links leading to a given spin 
have an identical sign. This leads to formation of magnetized 
clusters. 

A general Fokker-Planck equation that leads to the asymp- 
totic solution exp(— V (W,j3)) at x — > °°, has the form: 



d0» 1 . 
-dx~ = 2 di 



dV 



(8.3) 



where we define i = {i,j}, and differential operators <9j 
d/dWj. Differentiating the potential V, one obtains: 



dV 2W U 



dWij JijP 



■ tanh (Ri)- tanh (Rj) 



(8.4) 



A range of stochastic equations for the link noises can be 
obtained, by choosing different forms of the new diffusion 
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matrix In particular, we note that one may expect that 
a diffusion matrix that couples sites together over a distance 
of order of the expected correlation length might be expected 
to give a particularly efficient algorithm, as it tends to flip clus- 
ters of spins all of which have a similar spin orientation. For 
simplicity, we do not investigate this here, as we are interested 
in demonstrating a technique, rather than finding the most ef- 
ficient implementation. 



1. Constant diffusion 

For example, the simplest diagonal choice of 

%=/J/j% (8.5) 
leads to the following stochastic equation for the link noise: 



^ - - l -J-B-^ + t- (r) 



= -w ir 



JijB{mi + mj}+E,ij{T) , (8.6) 



where m t = tanh (/?,-) = tanh (/i/jS +Zj (Wy Q3) + Wjt (j3))), 
and: 

<gy (t) £ 7 (t')) = f}Jij8g8 jf 8 (t- t') . (8.7) 

Changing variables to R, = h{B + E/W^~Q3), with corre- 
sponding noises ^ — Y,j (§y + and an effective gain of 
gi = B £ j Jij, this reduces to 



dRj 



-Ri + gi tanh (Ri) + Bh, + £j3/y tanh (Rj) + £ (t) . 

(8.8) 

The important feature of this exact equation is that no addi- 
tional weighting is required. Each link noise equation is well 
localized, only scaling with the total lattice size. That is, for a 
D-dimensional lattice and nearest neighbour couplings, there 
are just MD link equations for M lattice points. The algo- 
rithm can be improved further by implementing link noises 
with variable correlation lengths for calculations near the crit- 
ical point, in order to spin-flip large clusters more quickly, and 
to reduce the problem of critical slowing-down. This could be 
achieved by having larger noise coefficients for longer wave- 
length Fourier coefficients. 

One can understand the equations physically as having a 
similar behaviour to the equation for the gain of a laser, with 
the first term causing loss and the second term gain, although 
with a nonlinear saturation as well. The first term is dominant 
at high temperature (small B), while the second term dom- 
inates at low temperature (large B). The external magnetic 
field term is like an injected field in the laser equations. The 
fourth describes correlations, while the last is a noise term. 



B. Example: 

As an example, consider the uniform two-node case again, 
where there is only one link and the two stochastic variables 




Figure 1: Stochastic calculation of two spin correlation: 7/3 = 1; 
4000 trajectories; step-size .05; semi-implicit method with 3 itera- 
tions. 




Figure 2: Stochastic calculation of two spin correlations over a range 
of J/3 ; comparison to exact results. 



are perfectly correlated. The stochastic equation is then: 
dR 

— = -R + 2j3/tanh (J?) + Bh + % (t) , (8.9) 
ax 



with 



(§(t)§(t')>=2 J 8/5(t-t / ) 
The correlation function is calculated from 



(afof> = ( [tanh (*)] 



(8.10) 



(8.11) 



The results of a simulation of Eq ( 18.9b are shown in Fig (Q~|). 
The corresponding correct result for the two-spin correlation 
is given by Eqs dT28i > and d733T > as: tanh(/j3) = tanh(l) = 
0.7616. Detailed results over a range of temperatures are corn- 
paid with exact results at thermal equilibrium in Fig (0. 

The sampling error in an ensemble of jV trajectories can be 
estimated as a j\fjY , where a is the standard deviation of the 
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Figure 3: Sampling error of two spin correlations. 
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Figure 4: Nearest-neighbor correlations for a 10 x 10 lattice as a 
function of inverse temperature. Circles: Results from stochastic 
calculations, 1000 trajectories. Solid line: Exact solution in the limit 
of an infinite lattice. 



calculated results, and assuming a nearly normal distribution. 
The actual sampling error for this simulation varies in time, 
and was estimated as 0.005, for large times - near equilibrium 
- as shown in Fig ©. 

Given this estimated error, the calculated stochastic result 
for the correlation agrees with the exact solution within the 



sampling error. 



C. Two-dimensional lattice calculation 

As a non-trivial example calculation, we consider a 10 x 10 
Ising model with periodic boundary conditions. Couplings 
are nearest neighbor, on a rectangular lattice with Jjj = 1 and 
h = 0. 

The numerically calculated nearest-neighbour correlation 
function is given for six different inverse temperatures j3. 
Once the relevant stochastic averages have reached steady- 
state, they are time-averaged as well as stochastically- 
averaged to give the correlation functions. 

The results are shown in Figure (3), along with a compar- 
ison to the known exact solution lBOll in the limit of an infi- 
nite lattice. The critical inverse temperature in this case is 
/3 C w 0.44, as seen in the exact solution. 



IX. SUMMARY 

We have shown how to obtain a general phase-space repre- 
sentation with positive-definite diffusion, for multiple SU(2) 
and more general SU(n) quantum systems, with couplings ob- 
tained from the corresponding operator algebra. In the case 
of qubits or two-level systems, the appropriate operator alge- 
bra is the spin half SU(2) algebra. This allows some further 
simplifications in obtaining evolution equations. 

The main application of these methods is to obtain stochas- 
tic methods for calculating either canonical ensembles or 
time-evolution of coupled atomic or spin systems. We have 
taken the exactly soluble Ising model as an example. The re- 
sulting stochastic equations were solved for correlation func- 
tions at finite temperature, and we found excellent agreement 
with known exact results. These techniques can also be ap- 
plied to more complex n-level cases, with time-evolution and 
coupling to external reservoirs. 
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